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Comparison  of  Numerical  Solutions  of  the  Vlasov  Equation 
with  Particle  Simulations  of  Collisionless  Plasmas 
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Washington,  D.  C.  20390 

and 

W.  L.  Kruer 
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Princeton,  New  Jersey  08540 

ABSTRACT 

This  paper  presents  numerical  solutions  of  the  Vlasov 
and  Poisson  equations  for  several  physically  significant 
problems  and  compares  these  solutions  with  the  results  of 
particle  simulations.  The  numerical  solutions  of  the  Vlasov 
equation  are  based  on  the  Fourier- Fourier  transform  method. 
The  spatial  representation  includes  up  to  85  modes  and  is 
capable  of  representing  strong  nonlinear  effects.  The  particle 
simulations  are  based  on  a  multipole  expansion  of  finite -size 
particles  about  their  nearest  grid  point  location.  Special 
techniques  are  used  to  suppress  the  noise  and  to  accurately 
control  the  initial  conditions  of  the  plasmas  so  that  quantita¬ 
tive  comparisons  with  Vlasov  solutions  can  be  made.  Close 
quantitative  agreement  between  the  results  of  the  two  simula¬ 
tion  techniques  is  observed.  The  problems  considered  are 


i 


one -dimensional,  with  periodic  boundary  conditions,  and  involve 
(l)  two -stream  instabilities,  with  equal  and  unequal  electron 
beams,  and  (2)  large -amplitude  electron  oscillations,  with  side¬ 
band  instabilities. 


PROBLEM  STATUS 

This  is  an  interim  report  on  a  continuing  problem. 


AUTHORIZATION 


NRL  Problem  H02-27 
DASA  HC-O^O 


I.  INTRODUCTION 


Numerical  simulation  of  collisionless  plasmas  may  be  achieved  by  either 
numerically  solving  the  Vlasov  and  Poisson  equations,  or  by  computing  the 
motions  of  a  large  number  of  charged  particles  that  are  moving  in  their  self- 
consistent  electric  field.  Although  the  two  methods  seek  the  same  ends,  they 
differ  fundamentally  in  their  approach  and  comparison  of  their  results  is  not 
trivial. 

This  paper  presents  numerical  solutions  of  the  Vlasov  and  Poisson  equa¬ 
tions  for  several  physically  significant  problems,  and  compares  those 
solutions  with  the  results  of  particle  simulations.  Close  quantitative  agree¬ 
ment  is  found.  Such  comparisons  provide  insight  into  the  validity  and  limita¬ 
tions  of  both  methods.  The  problems  considered  are  one -dimensional  with 
periodic  boundary  conditions,  and  involve  only  electrons  moving  over  a  uniform 
positively  charged  background. 

Numerical  solutions  of  the  Vlasov  and  Poisson  equations  have  been 
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carried  out  by  following  the  distribution  function  directly  in  phase  plane, 

and  by  transform  methods.  The  most  significant  transform  methods  to  date 

4 

are,  the  Fourier -Fourier  method,  in  which  the  distribution  function  is 

Fourier-transformed  with  respect  to  both  position  and  velocity,  and  the 

5 

Fourier-Hermite  method,  in  which  the  distribution  function  is  Fourier-trans¬ 
formed  v/ith  respect  to  position  and  its  velocity  dependence  is  represented  by 
a  series  using  Hermite  polynomials.  The  numerical  solutions  presented  here 
are  based  on  the  Fourier-Fourier  method.  Earlier  solutions  based  on  this 
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method  had  boon  limited  to  two  or  three  modes  in  the  spatial  representation. 
The  present  solo'-  mis  include  up  to  8**  modes  and  are  capable  of  representing 
strong  nonlinear  effects. 

Particle  simulation  of  plasmas  has  been  applied  to  a  variety  of  problems. 

In  this  method,  quantities  such  as  the  electric  field,  or  mean  velocity,  which 

depend  on  moments  of  the  distribution  function,  are  subject  to  random  noise. 

Special  techniques  are  used  in  the  present  solutions  to  suppress  this  noise  and 

to  accurately  control  the  initial  conditions  of  the  plasma  so  that  quantitative 

comparisons  with  Vlasov  solutions  can  be  made. 

Section  II  of  the  paper  reviews  the  Fourier -Fourier  method  used  in  the 

solution  of  the  Vlasov  equation,  which  in  its  main  features  follows  the  method 
-4 

ot  Knorr,  A  similar  review  of  the  particle  simulation  method  and  initializa¬ 
tion  techniques  is  given  in  Sec.  HI.  This  is  followed  in  Sec.  IV  by  the  results 
of  comparative  studies  of  the  two  methods  for  four  cases.  These  cases  were 
chosen  among  problems  which  had  been  considered  earlier  in  the  literature, 
so  that  comparison  could  be  made  not  only  between  the  present  solutions,  but 
also  with  earlier  numerical  studies. 

For  both  the  Vlasov  and  particle  solutions  presented  here,  time  is  meas¬ 
ured  in  units  of  u)^  ,  where  is  the  plasma  frequency,  length  is  measured 
in  units  of  the  periodicity  length  L  of  the  system,  and  velocity  is  measured 

in  units  of  Lto  .  It  follows  that  the  electric  field  is  measured  in  units  of 
P 

niLu)  /e .  where  e  and  m  are  the  electron  charge  and  mass,  respectively. 
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II.  FOURIER -FOURIER  TRANSFORM  METHOD 


In  terms  of  the  above  units,  the  one -dimensional  Vlasov  equation  for 
electrons  takes  the  form 


»f  9f  af 

9t  f  V  9x  '  E  9v  ' 


(1) 


where  f(x,v,t)  denotes  the  one -dimensional  electron  distribution  function  and 
E(x,t)  is  the  electric  field.  I^et  E(x,i)=  E  Xt(x,t)  +  Eint(x,t).  where  EtX* 
is  an  external  field  and  Emt  is  the  internal  field  due  to  electrons  and  the 
positively  charged  background.  The  internal  field  is  determined  by  Poisson's 


equation 


9E 


int 


9x 


=  l 


c 


f  dv 


Taking  Fourier  transforms  of  the  distribution  function  with  respect  to 
position  and  velocity  and  applying  periodic  boundary  conditions  in  space  with 
periodicity  length  L  =  l  yields 

l*"  r1 

H  (q,t )  =  \  exp(iqv)dv\  f(x,  v,  t)  exp  (-2zrinx)  dx 
_oo  o 

o 

The  functions  H  (q.t)  are  the  characteristic  functions  for  each  mode,  n 
n 

denotes  the  mode  number  in  space,  and  q  is  the  velocity  transform  variable. 
A  Fourier  transform  in  space  of  E(x,l)  yields  the  modes  of  the  electric  field. 


‘  f 


E(x,t)  exp(-2JTinx)  dx 


Since  f(x,v,t)  and  E(x,t)  must  be  real  valued,  we  have 
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H  (q,t)  =  H  (-  q.t)  , 
-n  n 


E  (t)  =  E  (t) 
-n  n 


After  transformation  and  truncation  at  a  finite  number  of  modes  m  ,  which 

max 

will  be  retained  in  numerical  computations,  the  Vlasov  equation  yields 


n  n  ~  , 

aT  +  Un  w  '■  2?cn,q't) 


where 


C  (  ,t)  =  -2jri  y 


>  E  (t)  H  (q.t) 

/  /  m  n-m 


is  a  convolution  term  which  comes  from  the  nonlinear  term  of  the  Vlasov 
equation.  Poisson's  equation  gives 


!  3T  i  Eint  =  -  P 


=  =  n  H  =  0) 

n  n  n 


for  n  ^  0 


and  Eint  =  0  . 
o 

Equations  (4)  are  solved  by  integration  along  their  characteristics,  which 

are  straight  lines  of  slop''  ZlTn  in  the  (t,  q)  plane,  as  shown  in  Fig.  1.  At 

each  time  step,  the  value  of  H  (q  ,  t)  is  obtained  from  the  iterative  formula 

n 

H^"+  \  q,  t)  =  H  (q-  ?./7nAt,  t  -  At)  +  ~ (q -  27Tn  At)  At  C  (q  -  27Tn  At,  t  -  At) 
n  n  47T  n 

+  ~  qAtC(l)(q  ,t)  ,  (6) 

47T  n  n 

in  which  the  superscript  denotes  the  number  of  iterations  carried  out.  The 


results  presented  in  this  paper  were  obtained  using  a  single  iteration. 
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From  the  definition  of  the  convolutions  C^,  Eq.  (6)  gives  the  improved 

approximations  for  n  =  m  m  ,in  terms  of  linear  com- 

n  max  max 

binations  of  the  preceding  approximations  H^.  The  matrix  01  the  coefficients 

n 

of  these  linear  combinations  is 


and  the  iterative  process  defined  by  Eq.  (6)  converges  only  if  all  eigenvalues 

X.,  for  j  =  -m  . +  m  ,  of  this  matrix  are  smaller  than  unity. 

j  J  max  max  1 

T 

Letting  A  denote  the  transposed  matrix  of  A,  we  have 

2  T  2 

|\j|  <  Trace  (A  A  )  =  (qAt)  mmax  u  . 

where  U  is  the  total  electrostatic  energy.  The  iterative  process  is  convergent 
therefore,  if 

WA'  (mmaxU)‘A  '  1  ■  <7> 

where  q  is  the  maximum  value  of  q  retained  in  the  computation.  This 
max  r 

criterion  is  satisfied  by  adjusting  the  time  step  At  according  to  the  magnitude 
of  the  electrostatic  energy. 

Values  of  q  in  Eq.  (6)  are  chosen  to  fail  at  grid  points,  as  shown  in 
Fig.  1.  The  values  of  Hn  at  q  -  27TnAt>  then,  do  not  fall  at  the  grid  points, 
and  must  be  interpolated  using  neighboring  points.  A  nine -point  Lagrangian 
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interpolation  was  used  in  the  computations  presented.  If  Aq  denotes  the 
interval  between  grid  points,  we  expect  the  distribution  function  f(x,v,t)  to 


be  adequately  defined  over  the  interval  -v  <  v  <  v  in  which 
'  max  max 

v  ~  l/Aq. 
max  '  M 

The  characteristic  functions  H  (q,t),  for  n  =  0 . m  are  evaluated 

n  M  max 

over  the  interval  -q  <  q  <  q  .  Values  of  H  (q,t)  for  negative  values 

max  max  n 

oi  n  are  found  from  the  reality  conditions  (3).  At  the  lower  boundary, 
q  -  -q  ,  the  values  of  H  (q  -  2ffn  At,  t  -  At)  needed  in  Eq.  (6)  are  unknown. 
These  values  are  set  equal  to  zero,  thus  introducing  a  cutoff  for  the  charac¬ 
teristic  functions  at  q  =  +q  .  The  introduction  of  cutoffs  m  and  q 

^  \nax  max  max 

is  equivalent  to  a  smoothing  of  the  distribution  function  f(x,v,t)  defined  by 


w  (v')dv'  \  w  (x')f(x  +  x',  v+v',t)dx' 
„  v  J  X 

-00  O 


where  f(x,v,t)  is  the  smoothed  distribution  function,  and 


w  (x)  =  i + 2  y 

x  L, 


cos  ZlTxvc 


(v) 


sin  q  v 
max 


are  weight  functions.  The  function  w  has  a  half -width  Ax  —  l/2m  ,  and 
6  x  max 

the  function  has  a  half  width  Av  —  The  choice  of  the  cutoff 

values  m  and  q  must  be  made  so  that  the  half -widths  Ax  and  Av 
max  max 

are  small  compared  to  characteristic  lengths  and  velocities  in  the  plasma. 


The  particle  code  considered  in  this  paper  uses  finite -size  particles 
having  a  Gaussian  charge  distribution  with  half -width  a  (see  Sec.  II1).  After 
Fourier  transformation  in  space  and  application  of  Poisson's  equa  ion,  this 
is  equivalent  to  multiplying  the  electric  field  of  each  mode  by  tne  form 

r  2. 

factor  W  =  exp[-(2iram)  J.  To  represent  the  effect  of  finite -size  particles 
in  the  Vlasov  solution,  this  form  factor  may  be  incorporated  into  the  convolu¬ 
tion  term: 


C  (q.t)  =  -2tri 
n 


w 

m 


E(t)  H  (q.t) 

m  n-m 


(8) 


m=-m 

max 


This  modified  form  of  the  convolution  reduces  to  the  original  form  given  by 


Eq.  (5)  if  the  form  factor  W  is  set  equal  to  unity. 

m 


The  introduction  of  the  form  factor  W  modifies  the  dispersion  relation 

m 


of  the  plasma,  which  becomes 


r<»  3  f  /3v 

- 2_ -  dv  = 

0)  -  27Tn  v 

-00 


27T 


0 


where  n  is  the  mode  number  and  U)  the  corresponding  complex  frequency. 
The  kinetic  energy  is  given  by 


(9) 


R 

where  HQ(q,t)  is  the  real  part  of  HQ(q,t)  and  the  electrostatic  energy  is 


given  by 


m 


max 


U 


■  l 


W 


n 


n=l 


(277n) 


H  (q  =  O.t) 


n 


(10) 
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terms. 


Each  term  of  the  convolution  array  C  is  the  sum  of  m 

n  max 

Since  there  are  m  terms  in  the  array,  the  computing  time  required  tc 

mAX 

evaluate  the  convolution  by  direct  summation  is  pioportional  to  .  The 

computing  time  is  significantly  reduced  by  using  discrete  Fourier  transforms 

to  evaluate  the  convolution.  The  arrays  H  (q,t)  and  E  (t)  are  Fourier- 

n  n 

transformed,  their  transforms  are  multiplied,  and  an  inverse  transform  of 

the  product  is  carried  out  to  obtain  the  convolution  array  C  .  Since  H  and 

n  n 

E  are  not  defined  cyclically,  but  instead  are  zero  for  n  >  m  ,  it  is 
n  7  1  rnax 

necessary  to  append  zeros  to  both  ends  of  the  arrays  before  carrying  out  the 
transforms.  Using  a  fast  Fourier  transform  algorithm,  ^  the  computing 


time  becomes  proportional  to  m  log,  m  .  This  method  is  advantageous 

max  2  max 

for  m  >10. 

max 


III.  PARTICLE  SIMULATION  METHOD 


The  particle  cods  investigates  the  electrostatic  behavior  of  a  plasma  by 
simply  following  the  motion  of  a  large  number  of  electrons  in  their  self- 
consistent  (and  any  externally  imposed)  fields.  Such  codes  are  widely  used  in 
computational  plasma  physics.^  ^  The  basic  scheme  consists  of  a  very  simple 
cycle.  The  positions  of  the  particles  determine  the  charge  density,  which  by 
Poisson's  equation  gives  the  self-consistent  electric  field.  In  one  dimension 
and  in  terms  of  the  dimensionless  units  used  in  this  paper,  Poisson's  equation 
takes  the  form 


i)E 


int 


N 


Ox 


■irl 


p.(x) 


j=l 
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where  p.(x)  is  the  charge  density  contributed  by  particle  j  and  N  is  the 
total  number  of  particles.  The  particle  velocities  and  positions  are  then 
updated  by  the  laws  of  dynamics,  using  a  standard  leapfrog  scheme 


v  = 


x  =  v 


Continuing  about  this  basic  cycle  advances  the  system  in  time.  Of  course, 

one  must  use  a  time  step  sufficiently  small  to  accurately  follow  the  time 

variation  of  the  forces  in  the  system.  One  tenth  of  is  generally  adequate. 

In  these  codes  one  also  inevitably  needs  to  discretize  space,  i.e.  ,  to 

introduce  a  regularly  spaced  grid.  How  one  then  defines  the  relevant  physical 

quantities  on  that  grid  is  the  principal  place  where  various  particle  codes 

12  13 

differ.  We  here  use  a  multipole  expansion  scheme  ’  and  finite-size  particles. 
The  charge  density  is  then  defined  on  the  grid  by  a  multipole  expansion  of  the 
particle's  charge  density  about  its  nearest  grid-point  location.  We  briefly 
illustrate  the  procedure:  Consider  a  particle  j  with  a  Gaussian  charge  distribu¬ 
tion  having  a  half -width  a. 


P.M  =  - 


exp[-(x  -x.)2/2a2] 
<2!T)1/2  a 


Here  x.  is  the  center-of-mass  location  of  the  particle.  Introduce  a  grid  and 

describe  the  particle  position  as  x.  =  n6  +  A*.  ,  where  6  is  the  cell  size  and 

J  J 

n  denotes  the  nearest  grid -point  location.  A*j  is  at  most  6/2.  Hence, 
assuming  6/2a  «  1,  we  expand  the  charge  density  as  follows: 


9 


P^X) 


expf-(x-n6)^/2a^]  ,  (x  -  n6) 

1/2  1  +  2 

(2  IT)  '  a  a 


Ax.  + 


Clearly  we  arc  replacing  the  finite -size  particle  centered  somewhere  in  the 
cell,  by  a  finite -size  particle  centered  at  the  nearest  grid  point,  plus  a 
finite- size  dipole  there,  plus  (in  principle)  higher-order  multipole  terms. 

In  practice  we  stop  at  the  dipole  correction.  Summing  over  a  collection  of 
particles  and  introducing  a  Fourier  transform  gives  the  total  charge  density 
in  Fourier  transform  space  as 


p(k)  -  -  exp  (-k^a^/2)  \  exp  (-ikn6)  [Q(n)  -  ikP(n)  +  .  .  .  ] 

n 

Here  Q(n)  and  P(n)  are  arrays  giving  the  net  monopole  and  dipole  moments 

2  2, 

associated  with  the  nth  grid  point.  Notice  the  form  factor  exp  (-k  a  /2), 
which  arises  due  to  the  finite  particle  size.  The  electric  field  is  now  determined 
simply  by  an  inverse  Fourier  transform.  The  force  on  the  particle  is  given  by 
the  same  multipole  expansion  procedure.  Physically  this  amounts  to  repre¬ 
senting  the  force  on  the  finite -size  particle  as  its  monopole  moment  times  the 
electric  field,  plus  the  dipole  moment  times  the  derivative  of  the  field. 

The  multipole  expansion  scheme  is  very  appealing.  First,  it  represents 

a  systematic  and  physical  way  to  introduce  the  spatial  grid.  Indeed,  an 

7  8 

expansion  parameter  has  been  exhibited.  Charge -sharing  schemes  ’  can  be 
related  to  the  multipole  expansion  by  stopping  at  the  dipole  approximation  and 
representing  derivative  terms  by  a  difference  over  cells.  Second,  the  multi - 
pole  expansion  scheme  relates  the  numerical  approximation  resulting  from 


10 


introducing  a  grid  to  physical  concepts.  One  is  investigating  the  physics  of 

a  plasma  of  finite -size  particles,  and  hence  there  are  some  modifications  of 

14 

the  plasma  behavior.  In  general,  the  finite  size  of  the  particle  enters  the 

analysis  via  the  form  factor,  the  Fourier  transform  of  the  particle  charge 

distribution.  For  example,  for  Gaussian  particles,  the  form  factor  is 
2  2  . 

exp  (-k  a  /2).  We  see  that  the  long -wavelength  (collective)  behavior  of  the 

system  is  essentially  unaltered,  but  the  short -wavelength  (ka  >  1)  behavior  is 

systematically  suppressed.  This  is  welcome,  since  short -wavelength  behavior 

(X  <  cell  size)  cannot  be  represented  accurately  due  to  the  finite  size  of  the 

grid.  Furthermore,  its  suppression  lowers  the  noise  level  and  hence  lowers 
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the  effective  collision  frequency.  This  yields  more  realistic  simulations  with 
fewer  particles. 

A  further  technique  used  for  reducing  the  noise  level  in  the  particle  code 
is  known  as  a  "quiet  start.  A  quiet  start  simply  refers  to  beginning  the 
calculat  on  with  ordered  initial  conditions.  Basically,  no  random  numbers 
are  used  to  set  up  the  calculations.  A  set  of  J  discrete  velocities  is  chosen 
using  the  probability  function 


P  =  P(v) 


■r 


f(v)  dv 


where  f(v)  is  the  desired  velocity  distribution  function.  The  procedure  is 
illustrated  in  Fig.  2  for  the  case  of  a  Maxwellian  distribution  function.  A 

set  of  J  values,  Pj  =  (j  -  0.  5)/J  with  j  =  1 . J,  are  chosen,  equally  spaced 

between  p  =  0  and  p  =  1.  The  corresponding  velocities,  v.  =  P  *(p.),  are  then 
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distributed  according  to  the  distribution  function  f(v).  Sets  with  J  varying 

from  100  to  1600  were  used  in  the  present  computations.  Equal  numbers  of 

particles  are  then  loaded  at  each  grid  point.  The  particles  may  be  loaded 

identically  at  all  grid  points  using  all  the  velocities  in  the  set  {  v }  ,  and  J 

particles  per  cell  are  then  needed.  However,  the  particles  form  a  number 

of  discrete  small  beams,  which  are  subject  to  instabilities,  and  spurious 

18 

oscillations  may  appear.  This  difficulty  is  reduced  by  using  a  larger  number 
of  beams.  To  achieve  this  without  increasing  the  total  number  of  particles, 
the  particle  velocities  are  staggered  so  that  each  discrete  velocity  is  repre¬ 
sented  only  at  every  kth  grid  point.  For  example,  the  set  of  velocities  may 
be  divided  into  k  subsets  by  picking  every  kth  value  of  the  original  set.  The 
particles  are  then  loaded  at  different  grid  points  with  different  subsets,  repeat¬ 
ing  the  same  loading  every  kth  cell.  The  repetition  length  k6  should  be  much 
smaller  than  any  physical  length  of  interest  in  the  problem. 

The  quiet  start  eliminates  the  initial  noise  level,  and  indeed  leaves  one 
the  option  of  starting  the  computations  with  specified  initial  conditions.  An 
initial  density  distribution  p(x)  may  be  obtained  by  giving  the  particles  initial 
displacements  from  their  uniform  distribution  at  the  grid  points.  If  £(x) 
denotes  the  displacement  of  a  particle  loaded  at  x,  and  p^  denotes  the  uniform 
density  before  displacements,  then  £(x)  is  found  by  integrating  the  equation 


dx  "  p(x  +  |) 


The  possibility  of  controlling  initial  conditions  is  important  for  detailed 

comparisons  with  other  codes— in  particular,  with  the  Vlasov  code  considered 
in  this  paper. 
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IV.  RESULTS 


Case  A 

Consider  a  two-stream  instability  resulting  from  the  initial  conditions 
defined  by  the  distribution  function 


f(x,  v,  t  =  0)  =  fQ(v)  [  1  +  2€  cos  27Tx] 


with 


f  (v)  =  - TTTT 

o  ,,.1/2  3 


(2ir)  '  v 


2  2,2 

v  exp  (-v  /2vth) 


(ID 


(12) 


th 


-2 


and  v  =  0.30 /it,  €  =  2.  5  X  10  .  These  initial  conditions  correspond  to  a 

system  length  L  =  10.  5  The  initially  excited  mode,  as  shown  in  Eq.  (11), 

has  a  wavelength  equal  to  the  length  of  the  system.  The  linear  growth  rates 

19 

for  this  problem  have  been  computed  by  Grant  and  Feix.  The  first  mode 

is  the  only  unstable  mode  and  has  a  growth  rate  y  -  0.25.  Vlasov  solutions 

20 

for  this  problem  have  been  carried  out  by  Armstrong  and  Montgomery 
using  the  Fourier -He rmite  method.  Comparisons  between  Fourier-Hermite 
solutions  and  particle -in -cell  solutions  have  also  been  made  by  Armstrong 
and  Nielson.  ^ 

The  solid  curve  in  Fig.  3  corresponds  to  the  total  electrostatic  energy  for 

the  Vlasov  solution,  with  m  =  21,  q  =  2  56,  and  Aq  =  4.  This  energy 

max  max 

grows  at  approximately  the  linear  growth  rate  from  t  =  10  to  t  =  20,  and 
saturates  at  2.2%  of  the  total  energy.  At  the  time  of  saturation,  the  bounce 
frequency  of  electrons  in  the  unstable  large  wave  (0.33)  is  roughly  equal  to 
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the  linear  growth  rate  of  the  instability  (0.14).  This  is  a  reasonable  result 
for  saturation,  since  a  large  wave  significantly  modifies  the  particle  dy¬ 
namics  in  a  time  -  l/ui  .  After  saturation,  the  electrostatic  energy  oscil- 

D 

lates  with  a  period  of  approximate!/  20.  The  frequency  of  trapped  electron 
oscillations  at  saturation  is  o>  =  0.31,  which  corresponds  to  a  trapping 

D 

period  Ttr=19.2. 

The  amplitude  of  the  first  mode  (n  =  1)  is  approximately  an  order  of 

magnitude  larger  than  the  amplitudes  of  the  other  modes  (n  ^  2).  The  higher 

modes,  however,  have  a  significant  effect  on  the  solution,  as  shown  by  the 

broken  line  in  Fig.  3  (which  corresponds  to  m  -  10).  The  Vlasov  solu- 

max 

tion  was  checked  by  reversing  it  at  t  =  20,  and  the  small  broken  line  near 

t  =  0  in  Fig.  3  shows  the  deviation. 

The  circles  in  Fig.  3  correspond  to  the  particle  solution  with  51,200 

particles  and  256  cells.  The  particle  half -width  is  a  =  1.256,  where  6  is 

the  cell  size.  The  initial  particle  positions  and  velocities  were  chosen  to 

match  the  initial  conditions  of  the  Vlasov  solution,  using  the  quiet -start 

method  of  Sec,  III  with  1600  velocities.  Since  there  were  200  particles  per 

cell,  the  initial  velocity  distribution  was  repeated  every  8  cells. 

Conservation  of  energy  was  checked  in  both  the  Vlasov  and  particle  solu- 

-4 

tions.  The  relative  energy  error  is  2  X  10  for  the  Vlasov  solution  and 

_4 

1.2  X  10  for  the  particle  solution.  To  give  some  idea  of  the  computing 

time,  the  unoptimized  Vlasov  code  required  -  13  min.  on  the  360/91  for  this 

problem.  This  compared  reasonably  with  the  time  of  -  15  min.  required  by 

22 

the  unoptimized  particle  code. 


14 


Several  additional  runs  on  the  particle  code  were  made  to  investigate 
spurious  oscillations  which  had  been  observed  in  earlier  runs  with  fewer 
discrete  beams.  The  solid  line  in  Fig.  4  shows  the  electrostatic  energy 
when  100  particles  are  loaded  identically  at  every  grid  point,  thus  forming 
100  discrete  beams.  In  contrast  to  the  smooth  behavior  of  the  1600-beam 
case  shown  in  Fig.  i,  strong  oscillations  of  the  fundamental  {n  =  I)  with  a 
frequency  CD  —  2  are  now'  superimposed.  These  oscillations  may  be  attri¬ 
buted  to  beaming  instabilities,  since  they  do  not  appear  when  the  number  of 
beams  is  increased.  With  the  present  quiet-start  technique,  the  most 
unstable  beams  are  those  representing  the  tails  of  the  distribution  function, 
since  they  are  the  most  widely  separated.  With  100  discrete  velocities,  the 

last  three  beams  are  located  at  v,  =  2.78,  2.99,  and  3.37  v  ,  respectively. 

b  tn 

The  frequency  of  beaming  instabilities  is  given  by  —  kv^,  where  k  is  the 

wave  number.  In  the  present  case,  k  corresponds  to  the  fundamental;  i.e.  , 

k  =  27T.  With  v  «  3  v.,  .  we  then  have  u>  —  1.  8,  which  agrees  with  the 
b  tn 

observed  frequency.  The  observed  growdh  rate  of  these  spurious  oscilla¬ 
tions  is  ~  kAV,  where  AV  is  typical  of  the  last  few  widely  spaced  beams. 

To  confirm  that  the  spurious  oscillations  are  caused  by  the  instability  of 

beams  around  v,  =  3  v  ,  ,  a  run  was  made  in  which  only  the  last  two  beams 
b  tn 

at  each  end  of  the  distribution  function  were  staggered.  The  resulting 
electrostatic  energy  is  shown  by  the  broken  line  in  Fig.  4.  Indeed,  the 
spurious  oscillations  now  appear  significantly  later  in  time.  This  is 
expected,  of  course,  since  the  less  separated  beams,  which  were  not 
staggered,  are  also  subject  to  instabilities  but  with  smaller  growth  rates. 
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Case  B 


We  now  consider  a  second  case  of  a  two -stream  instability  with  equal 
beams,  in  which  several  unstable  modes  are  allowed  to  grow.  The  initial 
conditions  for  this  case  are 

21 

f(x,  v, t  =  Oj  =  fQ(v)  [1  +  2€  cos  (Zltnx.  +  0^)}  ,  (13) 

n=l 

with 

f  (v)  =  - T75 -  [exp -{v+v  )2/v2  +  exp-(v-v  )2/v2]  ,  (14) 

°  2(jr)l/z  v  4  p  d  p 

P 

and  =  *^2  X  10  2,  v^-  2v^,  £  •  5  x  10~^.  The  Initial  phase  angles  are  chosen 

at  random  in  the  interval  (0,  2.7 r).  These  initial  conditions  correspond  to  a 

system  length  of  100  X^,  and  the  dispersion  relation  shows  that  modes 

23 

1  to  6  are  unstable.  This  case  was  studied  by  Morse  and  Nielson  using  the 
particle -in-cell  method. 

The  total  electrostatic  energy  for  the  Vlasov  solution  is  given  by  the  solid 

line  in  Fig.  5.  This  solution  was  carried  out  with  m  =  21,  q  =  1000  and 

max  max 

Aq  =  15.  6.  The  electrostatic  energy  reaches  approcimately  6.7 %  of  the  total 
energy.  We  again  observe  that  the  instability  saturates  when  the  electron 
bounce  frequency  in  the  dominant  wave  (0.  28)  is  approximately  equal  to  its 
linear  growth  rate  (0.26). 

The  broken  line  in  Fig.  5  corresponds  to  the  particle  solution  with  25,600 
particles  and  2  56  cells.  A  quiet  start  was  used  to  match  the  initial  conditions 
of  the  Vlasov  code,  including  the  same  values  of  the  phase  angles  0^.  Only 
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100  velocities  were  used.  The  discrete  beams  set  up  by  the  quiet-start  method 
are  rapidly  disrupted  by  the  collective  instability  in  the  present  case,  and  no 

spurious  instabilities  are  observed.  The  relative  energy  error  is 

-4  -3 

l.?.  X  10  for  the  Vlasov  solution,  and  0.8  X  10  for  the  particle  solution. 

Comparisons  of  densities  in  phase  space  for  the  Vlasov  and  particle  solu¬ 
tions  at  different  times  are  given  in  Figs.  6,7,  and  8.  For  the  Vlasov 
solutions,  numbers  from  1  to  9  denote  relative  densities.  Blanks  correspond 
to  densities  which  are  less  than  one  tenth  of  the  maximum  density.  Negative 
signs  correspond  to  negative  values  of  the  density,  which  occur  because  no  form 
factors  were  used  in  the  present  computation.  For  the  particle  solution,  an 
asterisk  was  printed  at  every  location  where  at  least  one  particle  is  present. 

The  results  of  Morse  and  Nielson  for  this  case  agree  qualitatively  with 
the  present  results.  The  electrostatic  energy  in  their  case  reaches  only 
5%  of  the  total  energy,  and  does  not  give  the  two  distinct  peaks  shown  in 
Fig.  5.  However,  the  case  considered  by  Morse  and  Nielson  corresponded 
to  a  longer  system  length  (L  =  500  X^),  In  addition,  these  authors  did  not 
control  the  initial  conditions  of  their  computations,  but  allowed  the  instability 
to  grow  from  the  noise  resulting  from  a  random  choice  of  initial  particle 
velocities.  Their  phase  density  plots  agree  with  those  of  Figs.  6,7,  and  8, 
showing  the  same  coalescing  of  the  eddies  when  their  system's  length  of 
500  X^  is  taken  into  account. 
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Case  C 


We  now  examine  an  instability  resulting  from  the  interaction  of  a  small 

beam  with  a  Maxwellian  plasma.  The  initial  conditions  are 

21 

f(x,  v,  t  =  0)  =  fQ(v)  [1  +  2e  ^  n  cos  (Zltnx  +  0^)]  ,  (15) 

n=l 

with 

yv>  =  — ~ — {np  exp(-v2/v^)  +  nb  exp  [-(v-vd)2/v2]}  .  (16) 

VP 

Here  v  -  10  2,  v ,  =  2. 6  v  ,  v,  =  0.25  v  ,  n  =0. 95,  n.  -  0.  05, 

p  d  p  b  p  p  b 

-4 

€  =  2.  5  10  and  the  initial  phase  angles  0^  are  chosen  at  random.  Thus 

the  small  beam  contains  5%  of  the  plasma,  and  its  mean  velocity  is  3.66 

thermal  velocities.  These  initial  conditions  correspond  to  a  system  length 

cf  100  X^.  The  dispersion  relation  shows  that  modes  1  to  9  are  now  unstable, 

23 

This  case  was  also  studied  by  Morse  and  Nielson,  using  the  particle -in¬ 
cell  method. 

The  total  electrostatic  energy  for  the  Vlasov  solution  without  form  factors 

is  shown  by  the  solid  line  in  Fig.  9.  This  solution  was  carried  out  with 

m  =  42,  q  =  2  5/v  ,  and  Aq  =  q  / 12 8 .  The  electrostatic  energy 

max  max  '  p  max' 

reaches  1.6%  of  the  total  energy.  Again  the  instability  saturates  when  the 
electron  bounce  frequency  in  the  dominant  mode  (0.16)  is  approximately 
equal  to  its  linear  growth  rate  (0.15). 

It  has  generally  been  found  that  each  trapping  region  requires  a  minimum 
of  8  to  10  modes  for  its  representation.  In  Case  A,  where  the  most  unstable 
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mode  is  the  fundamental  one  and  a  single  trapping  region  is  present  through 

the  system,  qualitatively  correct  results  were  obtained  with  m  =10. 

7  max 

In  the  present  case,  mode  m  =  5  is  the  most  unstable,  and  42  modes  were 

found  to  be  necessary  to  obtain  a  convergent  solution.  A  solution  with 

mmax  ~  ^  gave  qualitatively  different  results  beyond  saturation.  For 

CO  t  >  40,  the  m  =10  solution  gave  trapping  oscillations  which  continued 
p  max 

to  grow  in  amplitude  instead  of  dropping  to  the  low  levels  shown  in  Fig.  9. 

The  circles  in  Fig.  9  correspond  to  the  particle  solution  with  61,  440 
particles  and  2  56  cells.  A  quiet  start  was  used  again,  to  match  the  initial 
conditions  of  the  Vlasov  solutions.  In  the  present  case,  960  velocity  classes 

were  used  to  maintain  the  small-beam  instabilities  at  a  low  level.  The  rela- 

-4  -4 

tive  energy  error  was  —  1C  for  the  Vlasov  solution  and  —  3  10  for  the 

particle  solution. 

The  results  of  Morse  and  Nielson  in  this  case  again  agree  qualitatively 
with  the  present  results.  The  saturation  electrostatic  energy  in  their  solu¬ 
tion  was  2%  of  the  total  energy  instead  of  the  1.  6%  in  the  present  solutions. 
This  difference  may  again  be  due  to  the  longer  length  of  their  system  (200  A^) 
and  to  random  initial  perturbations. 

Case  D 

This  case  is  concerned  with  sideband  instabilities  resulting  from  the 
motion  of  trapped  particles  in  a  large -amplitude  electrostatic  wave.  The 
excitation  of  large -amplitude  waves  by  means  of  an  electrostatic  probe 
immersed  in  a  warm  plasma  has  been  described  by  Wharton,  Malmberg, 
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and  O'Neil.  These  experiments  showed  the  expected  amplitude  modulation 

of  the  wave,  which  is  attributed  to  electron  trapping,  but  they  also  disclosed 

the  appearance  of  sidebands  to  the  frequency  of  the  main  wave.  The  growth 

2  c 

of  these  sidebands  has  been  attributed  by  Kruer,  Dawson,  and  Sudan  to  an 

instability  due  to  particles  trapped  in  the  large -amplitude  wave,  and  has 

26 

been  observed  by  Kruer  and  Dawson  in  particle  simulations. 

In  the  present  computations,  the  initial  distribution  function  of  the  plasma 
is  defined  by 

42 

f (x ,  v ,  t  =  0)  =  f  (v)  [1  +  2c  /  ncos  (27Tnx  +0)1  ,  (17) 

o  i_t  rn 


n=l 


with 


£o<v)  3  ~T/Z  «P[-<v2/*4>] 

(2  7T )  ' 


(18) 


Here  v  =  1.06  /447T,  €  =  0.0002  and  the  initial  phase  angles  are  chosen  at 
random.  These  initial  conditions  correspond  to  a  plasma  length  L  =  130 


Mode  n  =  5  is  then  driven  from  t  =  0  to  t  =  6  by  the  external  field 
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0Xt 

E  (x,t)  =  E  sin  (oj  t  +  kx)  ,  (19) 

JJK  O 

with  E^„/v  .  =0.3  and  =  1.06.  The  driving  frequency  w  is  the  Bohm- 
DR'  th  o  °  o 

Gross  frequency  corresponding  to  mode  n  =  5.  and  the  ratio  of  the  phase 
velocity  of  the  driving  wave  to  the  thermal  velocity  is  ui^/Ztnv^  =  4.4. 

The  electrostatic  energies  of  the  main  wave  and  sidebands  from  the  Vlasov 
and  particle  solutions  are  compared  in  Figs.  10  and  11.  The  particle  solu¬ 
tion  was  carried  out  with  256  cells  and  200  particles  per  cell.  A  quiet  start 
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was  used  with  1600  discrete  velocities  to  represent  the  distribution  function. 

The  particles  were  given  displacements  to  match  the  initial  density  perturbation 

defined  by  Eq.  (17).  The  particle  half-width  was  a  =  26  =  A^,  where  6  is 

the  cell  length.  The  Vlasov  solution  was  carried  out  with  m  =  42,  thus 

max 

allowing  approximately  8  modes  for  each  trapping  region.  The  truncation  in 

velocity  transform  was  at  q  =  8/v  ,  and  the  grid  spacing  was  Aq  =  l/8v  ,  . 

max  th  th 

These  values  correspond  to  a  velocity  resolution  Av  —  ff/q  ^  —  v  /4  and 

a  maximum  velocity  v  —  l/Aq  =  8v  ,  .  A  form  factor  was  applied  to  the 

electric  field  of  the  Vlasov  solution  corresponding  to  the  particle  half -width 

-4 

a  =  A^  used  in  the  particle  code.  The  relative  energy  error  was  4  X  10 

-4 

for  the  particle  solution  and  3  X  10  for  the  Vlasov  solution. 

The  main  wave  energy  (electrostatic  energy  of  mode  n  =  5)  and  the  lower 
sideband  energy  (sum  of  the  electrostatic  energies  of  modes  n  =  1  to  4)  are 
shown  in  Fig.  10  on  a  logarithmic  scale.  The  main  wave  energy  rises  rapidly 
during  the  driving  period  to  29.4%  of  the  initial  kinetic  energy,  after  which 
it  oscillates  with  approximately  the  trapping  oscillation  period  (t  =  20.7). 
We  observe  close  agreement  between  the  vlasov  and  particle  solutions  for 
the  main  wave.  Indeed,  the  two  curves  are  nearly  identical  until  late  in  the 
simulation.  The  lower  sideband  energy  from  the  Vlasov  and  particle  solu¬ 
tions  grows  at  the  same  rate  (a  ten -folding  time  of  ~  31)  and  even  saturates 
at  the  same  level.  The  lower  sidebands  saturate  when  they  acquire  an  energy 
comparable  to  that  of  the  large  wave.  Then  the  sideband  waves  disrupt  the 
particle  trapping  in  the  original  large  wave,  as  confirmed  by  phase  space 
plots.  It  should  be  emphasized  that  this  problem  is  very  nonlinear  and  is  a 
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strong  test  of  both  of  the  simulation  techniques.  We  are  accurately  following 
not  only  sizeable  oscillations  in  the  large  wave  energy,  but  also  simultaneous 
growth  of  oscillations  at  other  wave  numbers. 

The  upper  sideband  energy  from  both  solutions  is  shown  in  Fig.  11.  The 
main  wave  energy  has  been  repeated  on  this  figure  to  provide  a  reference. 
Again,  the  two  solutions  agree  well  and  even  saturate  at  the  same  level. 

The  saturation  level  of  the  upper  sideband  is  approximately  an  order  of 
magnitude  below  the  saturation  level  of  the  lower  sidebands.  This  lower 
saturation  level  is  reasonable  since  the  upper  sidebands  have  phase  velocities 
less  than  that  of  the  main  wave  and  hence  are  more  readily  damped  by  the 
particles . 

Several  additional  solutions  were  carried  out  with  both  the  Vlasov  and 
the  particle  codes.  The  Vlasov  solutions  included  a  computation  with 
mmax  =  and  a  weaker  f°rm  factor  (a  =  \^/Z).  This  computation  showed 
only  minor  variations  from  the  results  of  Figs.  10  and  11.  The  particle 
solutions  included  one  with  fewer  beams,  to  describe  the  distribution  func¬ 
tion.  This  calculation  showed  that  details  of  the  saturation  levels  (but  not 
the  growth  rates)  are  sensitive  to  the  number  of  beams  used.  This  is  rea¬ 
sonable,  since  the  trapped  particles  responsible  for  the  sidebands  come 
from  the  tail  of  the  initial  distribution,  which  is  rather  poorly  represented 
if  too  few  beams  are  used. 
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SUMMARY 


This  paper  has  presented  quantitative  comparisons  of  particle  simu¬ 
lations  with  multiple-mode  solutions  of  the  Vlasov  equation  including  up  to 
85  modes.  Previous  solutions  of  this  type  had  been  limited  to  a  few  modes 
only.  The  problems  considered  ranged  in  complexity  from  a  two- stream 
instability  involving  a  single  unstable  mode  and  low  electrostatic  energy  (2.2% 
of  the  total  energy)  to  an  instability  due  to  particles  trapped  in  a  large-ampli¬ 
tude  plasma  wave.  By  using  quiet  starts  to  initialize  tha  particle  simulations 
and  using  a  sufficient  number  of  beams  to  suppress  beaming  instabilities, 
c’ose  agreement  was  found  between  the  two  methods. 

Since  the  two  methods  differ  fundamentally  in  their  approach,  the 
agreement  found  confirms  their  validity.  However,  the  problems  considered 
have  shown  limitations  in  both  methods,  which  must  be  taken  into  account  in 
the  physical  interpretation  of  numerical  simulation  results.  Discrete  particle 
effects  in  particle  simulations,  which  are  particularly  evident  in  regions  of 
low  density  in  phase  space,  yield  beaming  instabilities  which  must  be  mini¬ 
mized  or  accounted  for  in  the  physical  interpretation  of  the  results.  Similarly, 
solutions  of  the  Vlasov  equation  tend  to  develop  i.<  reasingly  fine  structures 
with  increasing  time.  The  fine  structures  are  suppressed  by  truncation  of  the 
Fourier  expansions  to  a  finite  number  of  modes,  but  enough  modes  must  be 

retained  to  make  the  half-widths  Ax  =  1/ 2  m  and  Av  =  7r/q  small 

max  max 

compared  to  the  characteristic  lengths  and  velocities  of  the  phenomena  being 
considered.  As  indicated  in  the  discussion  of  Case  III,  approximately  8  to  10 
modes  are  needed  to  represent  each  trapping  region  and  the  solution  may  he 
altered  in  its  general  character  if  fewer  modes  are  retained. 
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Fig.  10  -  Electrostatic  energy  of  the  main  wave  and  lower  sideband.  Case  D 


